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Starting from the total energy expressions within density functional theory, we are able to perform 
a comparison of several currently used charged-defect finite-size supercell correction schemes in a 
unified manner. This approach also provides a framework for a further development of corrections 
not only for DFT supercell calculations, but also for more advanced methods and for complex 
geometries. The comparison is performed for three separate defect cases: a gallium vacancy in 
GaAs, a beryllium interstitial in GaAs and a vacancy in diamond. We found two methods working 
sufficiently well for all three cases: a method which is very similar to one presented by Freysoldt, 
[l[ and a slightly altered potential alignment method. 
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I. INTRODUCTION 

The usual low concentration of defects in connection 
with the long-range Coulomb interaction is a difficult 
case for the supercell calculations within the density func- 
tional theory (DFT) framework. Especially if the defect 
is charged, the calculated formation energies strongly de- 
pend on the supercell size. Several correction schemes 
have been proposed to allow a calculation of defect prop- 
erties in an effectively low concentration by using only 
a small supercell. [ll, [3, H, 01 However, a consensus on 
the validity and applicability of each method seems to 
be missing, although numerical comparisons have been 
done. [1| 

In this paper, we start by deriving a generalized 
approach for the charged defect supercell calculations, 
which is based on the construction and comparison of 
the DFT total energy equations for the supercell and for 
a much larger cell. This is useful in several ways: a) 
This allows us to properly compare several contributions 
of the corrections and several previous correction meth- 
ods. We show how several schemes come out as limits 
or approximations of the general equations. We particu- 
larly concentrate on the method recently introduced by 
Freysoldt et al.,[l| and also the Makov-Payne and the 
potential alignment schemes Q are reviewed, b) It shows 
the required approximations in a clear manner. Even if 
the final results look intuitively simple, the way to get 
there has a few corners, c) Paves the way for the de- 
velopment of even better schemes. We present two more 
schemes in this paper, d) It should prove useful in de- 
veloping correction schemes for more advanced methods 
such as hybrid-functionals or GW , and for more com- 
plex geometries such as interfaces or clusters. The first 
three items are considered in detail in this article. The 



last one is briefly considered in the discussion part, but 
mostly left for later study. 

Finally, we compare these methods for three defects: 
gallium vacancy in GaAs, carbon vacancy in diamond 
and beryllium interstitial in GaAs. In order to find the 
results to compare at, we use a set of calculations with 
increasing supercell sizes to extrapolate the results to the 
low concentration limit. 



II. CORRECTION SCHEME 

The defect calculation is usually approached through 
the concept of the formation energy, defined as [sj] 
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-£;tot[bulk] - 

-q[EF + E, + AVl (1) 



where ii^tot l^"^] is the total energy of the supercell with 
the defect X (in the charge state q) and i?tot[bulk] is 
the total energy of the supercell of bare GaAs or GaAsN 
bulk, depending on the case. The chemical potentials 
fii of rii added or removed atoms allow us to describe 
various growth conditions. Here, Ey is the valence band 
maximum (VBM) at F-point in the bulk material, Ep 
is the Fermi energy with respect to Ey, and AV is the 
shift term used to align the potentials in between the two 
superceUs. 

Before moving on to the finite-size supercell considera- 
tions, we divide the comparison of charge defect and bulk 
cases in to two parts: comparison of a neutral defect to 
the bulk, and the comparison of a charged defect to a 
neutral defect 



Ef[X°] = E[X°] - S[bulk] - 



(2) 



'Electronic address: Ihannu.komsaQtut.fil 
t Electronic address: [tapio . rantala@tut . fll 



E^iX'^] E^[X°] + AE^iX'i] 
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where 

AE^iX"] = EiX"] - [E[X"] + qiE^X""] + Ep)] (4) 

We see, that m this approach the addition/removal of 
electrons, the VBM energy is in the energy reference of 
the neutral defect. The imphcations of this will be dis- 
cussed later. Moreover, Comparison of charge densities 
and potentials between charged defect and bulk can be 
difhcult, especially in the case of relaxed geometry. On 
the other hand, comparison of charged defect and neutral 
defect properties is often much easier. 

Naturally, defect calculations within DFT-LDA also 
suffer from the band-gap problem. This problem is im- 
portant to acknowledge in the analysis of the results, but 
does not affect the analysis of the supercell-size correc- 
tion methods. 



A. Charged-defect finite-size supercell corrections 

We will start by writing and comparing the total en- 
ergies for neutral defect, charged defect in a small super- 
cell, and a charged defect in a much larger supercell. For 
the neutral defect charge distribution po the total energy 
within the DFT framework is 

E[po] = T[po] + E,s[po] + E^cipo] (5) 

Adding a localized charge distribution 6q results in a re- 
distribution of the surrounding electrons. We circumvent 
this problem, by just writing the total density as po + Qd- 
Moreover, we denote the periodically repeating charge 
distribution of small supercell as qd- The total energies 
for the supercell of volume (l and a large cell of volume 
are then 

[po + Qd] = T'' [po + qd] + Eg [po + qd] + Eg [po + qa] 
E%o + qd] = T'^po + qd] + E'g[po + qd] + Eg[po + qd] 

Usually, one would like to calculate the formation energy 
of a single defect in an infinite crystal, or at least of much 
lower concentation of defects than is possible to obtain 
in the supercell calculations i.e., we would like to have 
E^^[po + qd]- Instead, what we get from the supercell 
calculations is E^[po + qd]- Thus, what we are looking 
for in here, is a correction AE, such that E^[po + qd] — 
E^^[po] required for Eq. [Hcan be obtained from E^^[X'^] — 
E^[X^] + AE. Also notice, in the following, we have 
taken the ionic configuration to be the same in all three 
cases. 

Due to the locality of the kinetic energy T and the 
exchange-correlation energy E^ci and the localization of 
qd in the supercell, we will now assume that 

[po + qd] = [po + qd] + [po] - [po] (6) 

and similarly for E^c, so it follows that T^[po + qd] — 
T^[po] = T^[po + qd] - T^[po]- Unfortunately, this can 



not be done for the electrostatic energy, but we can now 
write the total energy differences in the large and small 
supercell, and find, that their difference, which is the 
correction that we are looking for, depends only on the 
electrostatic energy difference 

= ^"[po + qd] - E''[po] - (^E''[po + qd] - ^""[po]) 
= Eg[po + qd] - Eg[po] - (^Eg[po + qd] - Eg[po]) 

The electrostatic potential corresponding to po is de- 
noted as Vo (consisting of the external potential and the 
Hartree potential). The electrostatic energy is (omitting 
dr^^ for brevity) 

Ecs[po]^ J PoVc^t + ^ J PoVh ^ J PoiVcxt + ^Vn) (7) 

where the ^ takes care of the double counting of electron- 
electron interactions. Moving on to the charged case, 
there is a change in the Hartree potential Vh VH+Vq/o 
and subsequently in the electrostatic potential Vq ^ Vg -I- 
Vqfo- The electrostatic energy is then 

E"[po + qd] = J^yo + qd + n){V,^t + ^Vh + ^Vq/^^) 

where n = —q/fl is the neutralizing background. V^/q 
is the change of electrostatic potential in small supercell, 
which is also, due to the linearity of Poisson equation, the 
solution for qd- In the periodic case, we have similarly 
added charge qd ~ qd, but per unit cell of volume f2, 
so the compensating background charge is n = —q/^, 
giving 

E'g[po + qd] - JApo + qd + n){V,^t + ^Vh + ^Vq0) 

Notice, that often h (and n) term is not included in DFT 
codes in the calculation of the electrostatic energy. How- 
ever, for the moment, we still leave this term in. 
It can be easily shown that 

f iqd + n)iV,^t + ^VH)= l{qd + n){V,^t + \vH) (10) 
and 

j qdVo - j PoVq,o (11) 

so that we are left with 

AE=^ f{qd + n)Vq/o -\j (qd + n)Vq/o (12) 

We have arrived at a rather intuitive form, which is 
often taken as a starting point for developing the defect 
correction formulae. Deriving this formula, qd was de- 
fined as the charge difference between the charged and 
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neutral defect calculations. In order to find more sim- 
ple correction formula, an analytic form is assumed for 
qd- In this case, however, one should also keep in mind 
that the corresponding potential needs to be correctly 
screened, which depends on the whole system. One can 
use the static dielectric constant and arrive at the Makov- 
Payne correction 0]. Alternatively, as was done in Ref. 
H) ^g/o can be divided in to the long- and short-range 
parts, where the screening only in the long-range poten- 
tial is handled with the static dielectric constant. 



1. Separation of long and short-range potentials 

In Ref. H, te long and short-range parts of the electro- 
static potentials were separated as 



9/0 



9/0 



(13) 



where V]^ is the potential solved from 



using the static 
dielectric constant. Similar separation is done for the 
Vg/Oi which is known from the calculation, and can then 
be used to determine the short-range potentials. Since 
qd is localized to the unit cell, Jj^ (idVqiQ can be changed 
to Jj^ IdVqiQ and a little rearrangement results in 

-AS = i ^(g<i+n)(V^,/o-"^^,/o)+^" 



^9/0 -^'^ 



v. 



9/0 

"(14) 

Which equals to Eqs. (6) and (7) in Ref. U, when n — > 0, 
except the second last term has an extra i. The minus- 
sign comes from the difference in the definition of the 
potential. Note also, that when given in this form, it is 
possible to calculate the correction corresponding to any 
concentration of defects (i.e., per volume fi). 

In calculations, the averages of potentials Vb, K?, and 
VqjQ are all set to zero, 
that it approaches zero as r 

T>sr _ 
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Next, we can choose Vq/q such 



Efly.'o+^as in Ref. 
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(X). Moreover, let's write 
With these consider- 



= o / ('Zrf + '5<z<i + n)(V;i'^-K"' + C) 



since 



/ K = / K -C = Jv.^o- V^' - C = / C (16) 



and h /j^ C = n Jfi^^j we get 
-AE ^ l 



{qd + Sqd + n)iV^'- -V^^) 
+ + + ^^^^ 

+ \n I < -\nl {Vl^). (18) 
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FIG. 1: The correction as a function of the width of the gaus- 
sian charge density. (Color online.) 



Finally, because /j^ fiV}^ = and \n ^ as 

— » oo, a very simple result is obtained: 
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• (qd + SqdW^^- -V^^-) (19) 
n 



We see that there is no dependence on the potential shift 
C, nor shifts A or B, but this is not really in disagreement 
with Ref. 1. In our approach, similar potential alignment 
term comes from taking the VBM with respect to the 
neutral defect energy reference (see Eqs. [2H4]). 



2. Makov-Payne 
Makov-Payne correction is 

2TTqQr 



iVn' + Vn ~ -n / (T/"- + yf)(15) where 



AE = q^al2erL 



(20) 



(21) 



is the second radial moment of the density difference. For 
SC lattice a = 2.8373 ^. When qd approaches delta- 
function, the Eq. [12] approaches Makov-Payne equation. 
This is demonstrated in Figure [1] for different supercell 
sizes and widths of qd- Similar calculations were also 
presented in Ref. 0. 

When qd is estimated as a gaussian distribution, the 
correction becomes somewhat smaller and the depen- 
dence on the supercell size decreases. Also note, that 
for a Gaussian qd, some of these terms can be calculated 
analytically (see e.g. Ref. |^.) 



4 



B. Potential alignment method 

In the potential alignment scheme, the difference of 
the electrostatic potential far from defect with respect 
to the bulk, AV, is read, essentialy resulting in energy 
correction AE = qAV. 

We calculated AV for a simple model charge density 
of a (nearly) point charge and a neutralizing background 
in a simple-cubic supercell geometry, expecting to obtain 
a simple dependence of AV on dielectric constant and 
supercell size. Indeed, we find the following form for the 
potential at the farthest point 

AV « 0.78-^ (22) 

where L is given in the units of bohr and the potential 
in the units of Hartree. The energy correction has again 
the same q^/SrL scaling although the constant is smaller 
than in the Makov-Payne correction for the cubic cell 
(2.8373/2 « 1.419). Later in this article, this energy 
correction is called analytic potential alignment. 

The similarity of the calculated and analytic potential 
alignments can be seen in Figures [U O and [H 



C. Energy comparison 

Here we consider how to best evaluate the VBM energy 
in Eq. [H In case there are no defect states near 
the VBM, there is no problem in the first place. However, 
even if there are defect states mixing with VBM, in the 
case of neutral defect, the electrostatic potential should 
converge to the bulk value fast (faster than r~^), and it 
can be obtained easily from the comparison of the po- 
tentials of the neutral defect calculation and of the bulk 
calculation. 

We find that taking the electrostatic potential differ- 
ence from the ion cores or the plane-averaged over the 
supercell far from defect gives very similar result. Espe- 
cially so in the unrelaxed geometries, but in the relaxed 
geometries the former probably proves easier to use. The 
usage of ^^^[X^] from the neutral defect calculation in- 
stead of bulk Ey is later denoted as VBM alignment. 
Notice the difference to the potential alignment method, 
where the potential difference is taken for each charged 
supercell of interest. 



D. Complete method and discussion 

Here we outline two correction schemes that seemed to 
work well for the studied cases. 

To sum up, the final correction scheme is then the fol- 
lowing: 

1. Obtain VBM alignment {Ey [X"]) from the compar- 
ison of neutral defect and bulk. 



2. Find gaussian qd which gives close match to the 
resulting change in the electrostatic potential be- 
tween charged and neutral defects. (When far from 
the defect.) 

3. Apply Eq. [Tni to get the final formation energy, 
along with the Eqs. [2H11 

This is the scheme I. 

We also found, that the following scheme (scheme II) 
seems to work surprisingly well: 

1. Obtain the VBM ahgnment {Ey[X°]) as before 

2. Correct the formation energies by using the analyt- 
ical form for the potential alignment Eq. 

III. APPLICATIONS 

In all of the calculations, we use planewave density- 
functional theory code VASP within the PAW-LDA for- 
malism. 9, 10, 11] In GaAs calculations, we have chosen 
gallium to have 3d frozen in the core and 400 eV cutoff. 
In diamond calculations, the cutoff is 500 eV. We use 
64-atom supercells with 4x4x4 k-points, 216-atom su- 
percells with 2x2x2 k-points, and 512-atom supercells 
with 2x2x2 k-points. 

We will now consider three test cases: a gallium va- 
cancy in GaAs, a vacancy in diamond, and a beryl- 
lium interstitial in GaAs. All of these are calculated in 
both the unrelaxed and the relaxed geometries. In all of 
the following formation energy figures, lines of the form 
aL^^ + bL^^ + c are fitted to the calculated values to ob- 
tain extrapolated values in the limit L^^ {L —^ oo). 

A. Gallium vacancy in GaAs 

The application of several correction schemes for the 
gallium vacancy are shown in Figure [H Makov-Payne 
scheme seems to work well, potential alignment underes- 
timates the formation energies. Scheme I tends to over- 
estimate the formation energies about 100 meV in large 
supercells and somewhat more in the 64-atom supercell. 
With VBM alignment, formation energies converge to a 
lower value. Consequently, scheme II energies are slightly 
lower than the extrapolated value from the uncorrected 
energies. However, the energies have very little variation 
over the supercell sizes in this scheme. 

B. Carbon vacancy in diamond 

Similar to Ref. [l|, we calculate the carbon vacancy in 
diamond for neutral, +2 and —4 charge states. In this 
case, due to the large band gap of diamond, there is no 
energy overlap of the defect states and the band edge 
states. This, along with large variations in stable charge 
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FIG. 2: Corrections for Vq^ in unrelaxed (left subfigure) and relaxed (right subfigure) cases. Color coding in left figure: no 
corrections (blue, solid), VBM alignment (blue, dashed), Madelung correction (red, solid), Madelung with VBM alignment (red, 
dashed), analytic potential alignment (green, solid), and calcuted potential alignment (green, dashed). Color coding in right 
figure: no corrections (blue, solid), analytic potential alignment with VBM alignment (scheme II) (green, solid), calculated 
potential alignment with VBM alignment (green, dashed), gaussian charge correction with VBM alignment (scheme I) (red, 
dashed) gaussian charge correction with charged cell VBM alignment (magenta, dashed). (Color online.) 



states, makes the case particularly suitable for studying 
finite-size supercell interactions without having to worry 
about the band-gap errors. 

The formation energies are shown in figure [3l Espe- 
cially for the +2 case, Makov-Payne notably overcorrects. 
Potential alignment seem to work fairly well. Scheme 
II surprises again, with very little variation and ener- 
gies close to the extrapolated value. Note, however, that 
for the Vj^ defect the analytic and calculated potential 
alignments differ considerably. 

Scheme I overcorrects again, but more than in the 
case of gallium vacancies. A reason can be traced to 
the charge distribution difference among the charged 
and neutral cases. The charge distribution difference 
of the relaxed Vj^ and V^^ with respect to the bulk 
case are shown in Figure [HJ The charge difference seems 
good in the neutral case, but in the charged cases, the 
added/removed charge is almost completely delocalized, 
except that e.g. Vq^ does not converge to —4, but —3.5. 
As the Poisson equation solution is also governed by the 
charge distribution away from the defect (large volume 
at large distance r), it seems worth trying to use these 
charges in the correction schemes. This fixes nicely the 
formation energies produced by scheme I, and also im- 
proves Makov-Payne results, except for the unrelaxed 64- 
atom supercell result. 

In any case, these charge distribution graphs seem to 
contradict with the basic premise of the correction meth- 
ods that there is a localized charge. Still, the corrections 
work fairly well. The occupied/emptied state was local- 
ized around the defect (not a host band), meaning that 
there must be a compensating change in the valence band 
electron density. This behavior in the charge differences 
was found in all three defect cases. 



The potentials V^^' and Vj^ and the calculated potential 
are shown in Figure SI We see that the potentials match 
very well when compared to the neutral case. This justi- 
fies the division of formation energy calculation into two 
parts. 



C. Beryllium interstitial in GaAs 

In order to test the method with something other than 
a vacancy, we calculated beryllium interstitial in GaAs. 
It was chosen, because we know from our previous studies 
p^ . that the defect states are well localized with very 
little geometric distortion. The formation energies are 
shown in Figure [6] General features for the neutral and 
charged cases are similar to the gallium vacancy case. 

Makov-Payne corrections work even better than in the 
case of gallium vacancy. Potential alignments without 
VBM correction tend to somewhat underestimate forma- 
tion energies as before. Here, scheme I works really well, 
and scheme II also relatively well even if it extrapolates 
again to a lower value. 

It seems, that using only the calculated potential 
alignment already corrects the formation energies about 
halfway, giving a fair confidence on the results reported 
in our previous study of beryllium defects, 

Comparison of the relaxed and unrelaxed cases for all 
defect types (and as can be seen in Figures [2] and show 
no systematic difference in the applicability of the correc- 
tion methods. During the derivation, these effects were 
approximated to be small, and this would indeed seem 
to be the case, although these defect cases were chosen 
especially chosen to show no major relaxation effects. 
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FIG. 3: Corrections for unrelaxed Vj'^ (left subfigure) and unrelaxed Vq^ (right subfigure) cases. Color coding as in Figure 
[21 except for the dash-dotted lines: Makov-Payne with effective charges (magenta, dash-dot), scheme I with effective charges 
(red, dash-dot). (Color online.) 
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FIG. 4: The xy-averaged potentials from 512-atom super- 
cell, similar to Fig. 2 in Ref. [J, of Vq^ against bulk (left) 
and against Vc (right). Blue solid: from calculation. Red 
dash-dotted: analytic periodic. Green dashed: the difference. 
Magenta dotted: analytic aperiodic. (Color online.) 



FIG. 5: The xy-averaged charge difference (multiplied with 
cell volume) of Vj'^ and Vq'* to the bulk. Unrelaxed geom- 
etry (left) and relaxed geometry (right) with horizontal lines 
roughly denoting the charge difference far from defect. (Color 
online.) 



IV. DISCUSSION 

Another problem with the extrapolation method is the 
amplification of errors in the set of calculations. If one 
of the calculations has an error of 100 meV for any rea- 
son, it can result in 1 eV difference in the extrapolated 
formation energy. Thus, it is advantageous to get rid of 
the extrapolation scheme. Taking this idea even further, 
a working correction scheme might allow to use lower 
precision in the calculation. One possibility would be 
to use coarser k-point meshes, although naturally other 
properties might also degrade. For example, it has been 
reported that coarse k-point mesh can lead to incorrect 



geometries. [T1,[T1,[T3 

Extending these correction schemes to more complex 
geometries such as interfaces requires again a proper 
model for the screening. This is straightforward enough. 
Alternatively, one could try to divide the space even fur- 
ther and write the total energy as the sum of the contri- 
butions in these regions. So far, very limited number of 
studies has been performed on defects at interfaces and 
warrants further investigations. 

When developing a correction for e.g the GW method, 
to the first approximation, the same trick could be done 
for the self-energy as was done here for the XC-energy, 
eventually yielding the same corrections. Of course, GW 
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FIG. 6: Corrections for Hef^ in unrelaxed (left subfigure) and relaxed (right subfigure) cases. Color coding as in Figure [2] 
(Color online.) 



provides a proper dielectric function for the screening 
which could be taken advantage of. Unfortunately, as 
GW calculations usually do not provide total energies, a 
somewhat different approach is probably needed. 

V. CONCLUSIONS 

A comparison of several previously introduced defect 
correction methods and a few new ones are compared in 
both analytical and numerical levels. First, we divide 
the problem of comparing charged defects in to bulk into 
comparison of neutral defect to bulk and further compar- 
ison of charged defect to neutral case. Then, by explicitly 
writing the total energies for a neutral defect, a charged 
defect in a small supercell, and a charged defect in a 
large supercell, it is possible to inspect the approxima- 



tions hidden in each method. This framework should also 
prove helpful in future development of methods for more 
advanced methods and for more complex geometries. 

We found that the method introduced by Freysoldt [l| 
works generally well. Moreover, during the inspection 
of the potential alignment method, we found a method 
which also worked surprisingly well, even if formally looks 
as just a scaled Madelung-energy. 
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